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Abstract 

In this paper, an analogy between the mathematical 
modeling of transonic potential flow and the flow in a 
cavitating bearing is described. Based on the similarities, 
characteristics of the cavitated region and j ump conditions 
across the film reformation and rupture fronts are devel- 
oped using the method of weak solutions. The mathemati- 
cal analogy is extended by utilizing a few computational 
concepts of transonic flow to numerically model the 
cavitating bearing. Methods of shock fitting and shock 
capturing are discussed. Various procedures used in tran- 
sonic flow computations are adapted to bearing cavitation 
applications, for example, type differencing, grid trans- 
formation, an approximate factorization technique, and 
Newton’s iteration method. These concepts have proved 
to be successful and have vastly improved the efficiency 
of numerical modeling of cavitated bearings. 

Introduction 

Cavitation in fluid film bearings, though recognized as 
early as 1886 when Reynolds introduced the theory of 


lubrication, is still a subject which draws intense debate 
as to the nature and mechanism of the phenomena. Vari- 
ous theories and conditions for cavitation have been put 
forward. However, only the collective works of J akobsson 
and Floberg (1957) and Olsson (1965), now known as 
JFO theory, have provided insight into the subject, which 
is both consistent with mass conservation and the physics 
of the problem. When the boundary conditions developed 
in JFO theory are applied to the Reynolds equation, the 
extent of cavitated regions and the performance of bear- 
ings can be predicted more precisely than by any existing 
method. This theory has yielded results which are in good 
agreement with experimental data. 

The subject of gas dynamics gained immense research 
interest around the turn of this century. This effort helped 
promote the development of supersonic aircrafts and large 
thrust rocket nozzles. Application of gas dynamics prin- 
ciples include turbine flows, gas lasers, aerodynamic 
windows, missile aerodynamics Jet engines and the flow 
around a body entering the atmosphere (Emanuel, 1986). 
Mach number M, a nondimensional parameter, which is 
the ratio of the flow speed to the local speed of sound, is 
the indicative index as to whether the flow is subsonic 
(M < 1), sonic (M = 1), or supersonic (M > 1). It is also a 
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measureof the compressibility of the flow. Inaconverging- 
diverging duct, the flow can range from subsonic to 
supersonic or the reverse. The converging section is called 
a nozzle when the flow is subsonic and is accelerating; it 
is called a diffuser when the flow is supersonic and is 
decelerating. The diverging section is called a diffuser 
when the flow is subsonic and is decelerating and is called 
a nozzle when the flow is supersonic and is 
accelerating.When subsonic and supersonic flow regimes 
exist, the flow is called transonic. It is also possible for an 
internal flow to be totally subsonic or supersonic through 
the nozzle. However, for a transonic flow to exist, a duct 
with a throat is essential. 

Prior to 1965, computational methods were rarely used 
in aerodynamic analysis and importance was placed on 
expensive and time consuming wind tunnel experiments. 
With the emergence of powerful computers, computa- 
tional aerodynamics has greatly advanced to the extent 
that the flow pattern past entire aircraft in different flight 
regimes can be predicted (Jameson, 1987). Such rapid 
growth in computational techniques can be attributed to its 
direct application in the design of aircraft and space 
vehicles. In addition, there was little recourse aside from 
expensive experimentation, due to the nonlinear nature of 
the governing equations which made them intractable to 
analytical modeling. Developments in computational 
methods applied in the lubrication area have been 
comparatively slow to evolve. In fact, to date, the numeri- 
cal algorithms developed by Elrod (198 1) and Kumar and 
Booker (1989) are the only effective numerical tools 
available in the analysis of cavitation in bearings. 

Transonic flow theory and the theory of lubrication are 
two distinctly different fields as far as the physical phe- 
nomena are concerned. While transonic flow theory deals 
with compressible fluid flow near sonic velocity, classical 
lubrication theory is generally concerned with the flow of 
a highly viscous incompressible fluid with a Reynolds 
number that is very small. However, due to the existence 
of subsonic and supersonic flow regimes in a converging- 
diverging nozzle and the existence of full film and cavitated 
regions in a bearing with a converging- diverging wedge, 
there exists a striking resemblance in the mathematical 
modeling of these two problems. Such an analogy can 
substantially benefit either field by suitably incorporating 
the advancements from one field to the other. 

In this paper, a mathematical formulation for a cavitated 
bearing is derived and compared with that of transonic 
potential flow. An analogy between these formulations 
are developed. The analogy is utilized to employ the 
method of characteristics and the method of weak solu- 
tions from transonic flow theory to determine the 
characteristics of a cavitated region and the jump condi- 
tions that apply across both film reformation and rupture 
fronts. The method of determining the location of film 


reformation using shock fitting and shock reformation 
techniques are discussed. A brief discussion of several 
techniques that are widely used in current transonic flow 
computation is provided. These techniques have already 
been suitably modified and incorporated into the analysis 
of cavitation in bearings. 

Mathematical Modeling 


Cavitated Bearing 


The conservation of mass flow, within the clearance 
between the stationary and moving surfaces of a bearing 
can be written, by lumping across the film thickness, as 


dph 

lh 


+ V - m = 0 


( 1 ) 


In the converging wedge of the bearing, the film thickness 
diminishes and the pressure is developed. In this region, 
the mass flux ntf can be represented by 




( 2 ) 


The first term on the right side is the mass flow due to 
shear (Couette flow) and the second term is the flow due 
to pressure gradient (Poiseuille flow). Somewhere in the 
diverging wedge of the bearing, the film ruptures and a 
cavitation region is formed which continues until the film 
is reformed again. In this cavitated region, the pressure 
remains essentially constant at the cavitation pressure and 
the mass flows across this region due only to shear along 
the film stnations. in the cavitated region, the film occu- 
pies only a portion of the volume, the remaining portion 
being filled by air, gas, or vapor. Mass flow in this region 
is given by 



(3) 


where 0 is the partial film content in the cavitated region. 
Although the film consists of incompressible fluid, if it is 
assumed that the density of the film varies due to the 
applied pressure, then the variable 0 can be provided with 
a dual interpretation 

0 = Jp/Po in the full film region where 0 > 1 
|v f /V t , in the cavitated region where 0 < 1 

where Vf and V t are the volume occupied by the film and 
the total volume, respectively. In the full film region, the 
variation of density will be governed by the bulk modulus 
of the liquid, that is. 


»- 


*dp ae 


(4) 
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Equation (4) also enables one to represent the pressure in 
terms of the density (or 0), and in essence, acts as the 
equation of state of the lubricant 
Mass flow through the entire bearing can then be written 
as 


U U A 

m = — p c h0 - g 
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( 5 ) 


where g is a switch function, which is introduced to 
remove the flow due to pressure gradient within the 
cavitated region and is defined by 


g = 


1 when 0 > 1 
0 when 9 < 1 


For a finite bearing, the flow due to shear occurs only in the 
circumferential direction while the flow due to pressure 
gradient is present in both the circumferential and axial 
directions. Hence, the two dimensional form of equation (1 ) 
can be written as 

dp c 0h + _3_rpchU0 p c Ph 3 30 'I 

dt + dx { 2 12*i g 3x J 


dz 


PcPh 3 ae > 

12p 8 dz 


= 0 (6) 


Equation (6), which can be used to describe the mass flow 
through the entire bearing, was developed by Elrod and 
Adams (1974) and termed a ‘universal equation’. In the 
full film region, equation (6) may be written as 


dE U dE 

dt 2 dx 



where E = 0h and K= -Ph 3 g/12*i. Equation (7) is an 
elliptic partial differential equation. The form of equa- 
tion ( 6 ) in the cavitated region is given by 



U dE 

2 dx 


= 0 


( 8 ) 


and is a hyperbolic partial differential equation. This may 
be easily demonstrated by differentiating with respect to 
t to get 

d^E U _d^E_ = d^E _ Ul d^E _ (q , 

dt 2 2 dx dt 3t 2 4 dx 2 ~ 

Equation (9) isacanonical form of the wave equation. The 
characteristic form of this equation has two real roots, 
±U/2. 

In the full film region, pressure increases to a maximum 
value and then gradually decreases until the pressure and 
its derivative simultaneously vanish, at a location where 
the film ruptures. The air/gas/vapor strips in the cavitated 
region begin with a pointed shape. When the film reforms, 
depending upon the upstream conditions, the reformation 
occurs abruptly. This is because the cavitated region is not 
able to signal the impending conditions to the upstream 
flow. Figure 1 represents typical profiles of pressure and 
fractional film content for a submerged journal bearing, at 
a particular axial location. The abrupt changes in 0 and 
the pressure gradient at the reformation boundary and the 
gradual changes at the rupture boundary can be clearly 



Figure 1.— Typical pressure and fractional film content distributions 
in a journal bearing. 
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seen. When the effects of cavitation are not formally 
considered, it is usually assumed that the film ruptures 
at the minimum film thickness and reforms at the maxi- 
mum film thickness. When the cavitation boundaries are 
determined, it is found that the film extends slightly 
beyond the minimum film thickness into the diverging 
portion of the bearing and, depending upon the lubricant 
supply conditions, the film reformation can occur at or 
around the feed groove. 

Transonic Flow 


The flow of a compressible fluid in thermodynamic 
equilibrium is governed by the Navier-Stokes equations. 
For a two-dimensional flow these equations can be written 
in vector form as 


dw ^ 3f + 9g _ 3R 3S 

dt dx dy dx dy 


( 10 ) 


where w is the vector of dependent variables; density, 
Cartesian velocity components, and total energy; f and g 
are the convective flux vectors; and R and S are the 
viscous flux vectors. Because the full Navier-Stokes 
equations are quite complex , approximations are generally 
made. One such simplification consists of assuming no 
viscous dissipation and that flow is irrotational. Conse- 
quently, equation (10) can be written in a quasi-linear 
form in terms of the velocity potential <(>. 
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It should be noticed that, in a subsonic flow regime 
(u,v < c), the coefficients of the second order terms will be 
positive and, in a supersonic flow regime (u,v > c), the 
coefficients become negative. This variation results in the 
equation being of the elliptic type with two imaginary 
roots in subsonic regions and of the hyperbolic type with 
two real roots in supersonic regions. 

When a subsonic flow slows down, it does so gradually. 
On the other hand, a supersonic flow, which can also 
decelerate gradually, normally slows down abruptly. 
Because the fluid in a supersonic flow is unable to signal 
the upstream flow of any flow or geometric changes. This 
is a typical characteristic of phenomena governed by 
hyperbolic-type equations. The abrupt changes lead to 
discontinuities in the flow which are known as shock 
waves and are a distinct feature of a supersonic flow in 
establishing the overall nature of the flow field. Of course, 
an internal flow can also emerge as supersonic without the 
presence of any shock if the outlet conditions permit 
Since viscous effects are neglected in the potential equa- 
tion, for internal flow, sonic conditions occur at the throat 
of the duct. For real fluids, the sonic line extends slightly 
downstream of the throat into the diverging portion of the 
duct. 

Computational techniques for potential flow have been 
extensively developed, since, extremely inexpensive 
computation is achieved with this formulation. Moreover, 
shock capturing and convergence acceleration techniques 
developed for potential flow have been found to be trans- 
ferable to more complex models using Euler equations. 
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Figure 2. — Transonic flow and cavitated bearing. 
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Analogy 

The similarities between internal transonic flow and 
cavitated bearing modeling are evident from the previous 
two sections. Figure 2 illustrates these similarities. Both 
phenomena are governed by similar mathematical formu- 
lations in the different regions and both have an embedded 
hyperbolic region within the elliptic region. The subsonic 
portion of the flow (M < 1) is analogous to the full film 
region (0 > 1). The sonic line (M = 1) is analogous to film 
rupture (0 = 1). The supersonic portion of the flow (M > 1) 
is analogous to the full film region (0 < 1). Also, a shock 
and film reformation have similar characteristics. It should 
be recognized that the Mach number and the inverse of 
fractional film content have similarities. The sonic line/ 
the film rupture locations are strongly influenced by the 
geometry of the flow, while the shock wave/film reforma- 
tion locations are due to the upstream conditions. The flow 
can also be fully elliptic or hyperbolic in both cases, 
although a completely cavitated bearing has no physical 
significance. Similar to compressible flows involving 
shocks, determination of the film reformation boundary is 
a difficult task. 

Although, it is seen that transonic flow and flow in a 
cavitated bearing have similarities, it should also be noted 
that they also differ in several respects. The essential 
differences between these two models are the following: 

( 1) In transonic flow, the type of the equation is changed 
due to the change in the sign of the coefficient of the 
second order terms; in the case of a cavitated bearing, the 
second order terms are totally lost in the hyperbolic region 
resulting in the reduction of the order of the governing 
equation. This sometimes causes oscillations at the 
boundary locations. 

(2) The entire flow is compressible in a transonic flow; 
but, in a cavitated bearing, the full film region flow 
although really incompressible is taken to be compress- 
ible and the cavitated region flow is of the compressible 
type. 

(3) The unknowns in an irrotational transonic flow are 
density and the potential function which are dependent on 
each other. On the other hand, for a cavitated bearing the 
only unknown is density (or 0). 

(4) Row in both Cartesian directions can occur in 
transonic flow, although the resultant velocity can be 
made to align with one of the coordinate axes by Jameson *s 
rotated difference scheme (1974). In bearings, generally 
the only motion is a direct result of the journal rotation and 
the velocity vector is normally taken to coincide with a 
coordinate axis. 


Having pointed out the basic analogy between the two 
models, extension of method of characteristics to deter- 
mine the path of the disturbances in the cavitated region 
and the method of weak solutions to determine the jump 
conditions across the discontinuities can be developed. 
Also, determination of the location of film reformation 
using shock fitting and shock capturing methods will be 
discussed. Since, our intention is to develop correspond- 
ing expressions for a cavitated bearing based on the 
analogy, the development details for transonic flow are 
not presented in detail, here. Interested readers are refered 
to, for example, Anderson et al. (1984). 

Method of Characteristics 

Hyperbolic equations have certain lines (or surfaces) 
which indicate the zones of influence and zones of 
dependence. The information about the flow is signalled 
along these lines which are called characteristics. This 
property is used to determine the value of the variable at 
a particular location in a hyperbolic region from the 
known value of the variable at a downstream location. 
This method of solving hyperbolic equations is the Method 
of Characteristics (MOC). 

Supersonic flow. - Assuming the free stream is aligned 
with the x axis, equation (1 1) can be written as 

(l - Mi )♦« + ♦„ = 0 (12) 

where is the freestream Mach number. In order to 
determine the characteristic direction, equation (12) is 
written in terms of the Cartesian velocity components 
along an arbitrary smooth curve C, and the determinant of 
the coefficients is set to zero along the curve. This will 
result in the differential equations for the characteristics. 
For this case, 


dy 

dx 


+ ! 

B 





(13) 


For a constant B, the equations describing the character- 
istics are obtained by integrating equation (13). This 
results in the following form: 


% = x - By 
B = x + By 


(14) 
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where % and t| are coordinates along the characteristic 
line. 

Cavitated region . - Within the cavitated region, the 
governing hyperbolic equation is written as 


This eliminates the need for the solution to be differen- 
tiable across the discontinuity. The mathematical theory 
of weak solutions for hyperbolic equations is a relatively 
recent development and may be utilized to determine the 
jump conditions across a discontinuity in a flow. 


dE U dE 

at + 2 dx 


(15) 


Along a curve x = x(t) in the x - 1 plane, E = E(x). For a 
particular curve x = x c (t), let dE = 0. 


Thus, 


Transonic flow. - Although the steady state formula- 
tions exhibit elliptic and hyperbolic type equations at 
different parts of the flow field, addition of an unsteady 
term results in hyperbolic type of equation. Consider a 
one-dimensional, scalar, nonlinear, hyperbolic partial 
differential equation 


dE 


3E 

aT dt + 


0E 

ax 


dx = o 


(16) 


du 

at" 



( 20 ) 


Using equation (15), results in 


U 3E ^ 3E , aE ( dx. U 
— — dt+ — dx,. = — ■ — 


2 ax 

Therefore, 


dx 


dx { at 2 ) 


dt = 0 (17) 


3x c U 

at ~ 2 


•\ 


§ =0 

dt 


(18) 


J 


The solution to the first of equation (18) determines the 
equation of the characteristics and the second one reveals 
the parameter that is constant along the characteristics. 
They are 


\ = *c 



E = 0h = constant 


(19) 


Since h is known, the value of 0 at any point in the 
cavitated region can be obtained from a point with a 
known value of 0 by tracing backwards along the char- 
acteristic line, Olsson (1965) discussed this characteristics 
approach in his treatise on dynamically loaded bearings. 


Method of Weak Solutions 


A genuine solution of a hyperbolic differential equation 
is one in which the dependent variable is continuous but 
discontinuities in its derivatives may occur. Alternately, a 
weak solution is genuine except along a surface across 
which the dependent variable is discontinuous. Only the 
scalar or vector form of first order and hyperbolic second 
order partial differential equations possess weak solu- 
tions. Since the dependent variable is not continuous 
across the discontinuity, an integral formulation is used. 


where u and F(u) are unknown variables. This can be 
rewritten as 


du 

dt 


A 

+ A 37 


= 0 


( 21 ) 


where A = A(u) = dF/du is called the Jacobian. Now, if 
w(x,t) is an arbitrary test function which is continuous and 
has a continuous first derivative but vanishes on the 
boundary and outside of an arbitrary domain D inthe(x,t) 
plane, then 


+ f^) w(x > t)dxdt - ° < 22 > 

When both u and F are continuous and have continuous 
first derivatives, it can be shown that 


JJ D (uw, + Fw x )dxdt = 0 (23) 

Functions u(x,t) which satisfy equation (23) for all test 
functions w are called weak solutions of equation (21). If 
the domain D contains a moving curve x(x,t), across which 
u is discontinuous as shown in figure 3 , equation (23) may 
be integrated by parts, utilizing equation (22), to get 

f ([u]cosai + [F]cosct 2 )ds = 0 (24) 

JZ 

where [ ] denotes a jump across the discontinuity and 
cos ctj and cos are the direction cosines normal to the 
discontinuity T (x,t). Since the integrand must vanish for 
all w, the condition for u to be a weak solution of 
equation (21) will therefore be 


[ujcoscq + [F]cosa 2 = 0 (25) 

Equation (25) determines the jump in the value of u 
across the discontinuity. This analysis is also valid for a 
system of equations in which case, u and F are vectors. 
The jump conditions. 
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Figure 3. — An arbitrary domain with a discontinuity. 


[U]£ = [F] (26) 

at 

where dx/dt is the velocity of the discontiuity, are 
popularly known as the Rankine-Hugoniot equations. For 
example, if u and F are vectors defined as 


u 


IT 


p,pu,pe + p— 


pu,pu + p,u 


f 2 ^ 

U 

pe + p— 


“|T 


then the Rankine-Hugoniot equations can be written as 


u„[p] = M 

u d[p u ] - [pu 2 + p] 

U d [pe + pu 2 / 2 ] = [u(pe + (pu 2 /2j + p)] 

where Uj is the velocity of discontinuity and e is the 
internal energy. The first two equations are called 
mechanical jump conditions. 


Cavitated bearing . - Consider the one-dimensional 
version of equation (6), which can be written as 


3E dm 
3t + 3x 


0 


(28) 


where m is the mass flux as defined by equation (5). 
Equation (28) can be written in a form similar to that of 
equation (21), that is, 


where A = A (E) = dm/dE. To proceed in the same manner 
as that used for the transonic flow illustration, it can be 
shown that for E to be a weak solution of equation (28), 
the jump condition to be satisfied at any discontinuity is 

[E] cos a x + [m] cos a 2 = 0 

At any time tj, the velocity of propagation of the discon- 
tinuity is determined by the tangent of the curve as shown 
in figure 4. The direction cosines are 


dt 


COS (X] = 

(*♦ 

therefore, 


j — and cosot 2 

dOa 

** x J t =,j 


cos ct) _ dx 
cos a 2 dt 


dx 



dt)2 

dx J»=M 


(31) 


Hence, equation (30) can be written in terms of the 
primitive variables as 


K - ( 0 , Or]£ - ItHl - ( 0h ) R ] 


12*i 


h ’ 8<8) s) L ■ ( h38<e) f ), 


= 0 (32) 


where L and R are the conditions at the left and right sides 
of the discontinuity, respectively. In other words, 


3E . 3E 


= 0 


(29) 


wf- = M 

is an equivalent Rankine-Hugoniot equation. 


(33) 
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Figure 4. — Path of a discontinuity. 


First consider the discontinuity due to film reformation. 
If it is assumed that the Film thickness variation across the 
discontinuity is negligible and that 0 r= 1 since full Film 
is formed at the right side of the discontinuity, equa- 
tion (32) for these conditions can be written as 


^ - 'E ■ 


>■ - - £ 


HI 


= 0 (34) 


If the reformation front is not moving with respect to time 
(analogous to a stationary shock), then the conditions for 
film reformation will be 


>■. ■ o ■ 



(35) 


Since (30/3x) > 0, obviously, 0 ^ 1. 


On the other hand, if the reformation front is in 
motion, equation (34) can be written as follows 


■" . — ](0 L . i) . . 

2 dt/ L ’ I2\i v. d* 


(36) 


Since (dO/dx) > 0, the right side of equation (36) is less 
than or equal to zero. Also, since 0 L <1, there can be two 
conditions 


(i) dx / dt > U/2 

For this condition, the left side of equation (36) is greater 
than or equal to zero. Hence, the only condition that will 
satisfy the equality is that both sides must be zero, that is, 
(30/3x)r = 0 and 0l = 1. This is the classical Reynolds 
boundary condition which is generally applied to deter- 
mine film rupture. 


(ii) dx / dt < U/2 


Now, the left side of equation (36) is also less than or 
equal to zero. Hence, q, cannot be uniquely determined 
from equation (36); it must be obtained using the charac- 
teristics. Equation (36) is treated as a differential equation 
of the following form to determine the velocity of the 
front: 


dx = U JJ_ 1 f 2 09^1 
dt 2 ' 12*i (1 - 0 L )l dxJ R 


(37) 


Although in transonic flow, the sonic line is not treated 
as a discontinuity; in the case of a dynamically loaded, 
cavitated bearing, there could be a discontinuity of 8 at 
the film rupture, when dx/dt>U/2 and there will be a jump 
in its gradient The jump conditions for this case can be 
developed in the same manner as previously described. 
The resulting expressions will be 


(i) dx /dt £ U/2 (de/dx) L - 0 and 0 R - 1 


(ii) dx/dt > U/2 


dx 

dt 


U 

2 


P 1 

i2n (e R - 1) 



These jump conditions are exactly the same as those 
determined by Olsson (1965). However, in that reference 
the conditions were derived using a mass balance across a 
fluid volume containing the discontinuity. 


Computational Treatment of Discontinuity 


Transonic flow . - The numerical computation of 
supersonic flow is complicated due to the presence of 
shock waves, across which the dependent variables and 
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their derivatives may be discontinuous. Two types of 
numerical techniques have been developed to analyze 
such flow fields and are known as shock fitting and shock 
capturing techniques. 

Shock fating technique . - This technique attempts to 
locate any discontinuities and treats them as boundaries 
between the regions of the flow field where regular solu- 
tions are applicable. Shock fitting is achieved by satisfying 
the Rankine-Hugoniot equations across the discontinuities 
while simultaneously ensuring that the solution on the 
downstream side of any shock is compatible with the rest 
of the flow field. The movement of the shock wave is 
obtained as a part of the solution. The flow field down- 
stream of each shock can be determined from free stream 
conditions. If the upstream conditions, initial shock slope 
and velocity are known, the shock acceleration and post 
shock pressure can be determined by combining Rankine- 
Hugoniot equations with the compatible equation. This 
technique is most convenient for governing equations 
written in nonconservative form. 

Several approaches have been devised to fit shocks 
(Moretti, 1974). The flow field is either partitioned by 
aligning any shock waves with grid lines or the 
discontinuities are treated explicitly, but not as bound- 
aries, in the computation. 

Shock capturing technique . - Unlike the shock fitting 
technique, with this method, the discontinuities are pre- 
dicted as a part of the solution without the requirement of 
any special treatment. By casting Euler equations in 
conservation-law form, the weak solutions and jump 
conditions are built in. The conservative form of govern- 
ing equations and the discretization automatically allow 
prediction of the shock wave speed and the strength 
(Lax, 1954). 

Because of the simplicity in approach, this technique is 
most popular in the computation of flows with shocks. 
This technique also has several variants which include 
flux splitting and split coefficient matrix methods. Shock 
waves predicted with this technique can be smeared over 
several grid spaces and the application of surface bound- 
ary conditions can be difficult. Due to the wavelike nature 
of hyperbolic equations, boundary errors are propagated 
into the flow field which results in instability in the 
computation. In general, shock capturing techniques are 
applied to predict internal shocks while boundary shocks 
aie fitted. 

Cavitated bearing . - Similar to the previously described 
approaches for transonic flow, discontinuity fitting and 
capturing techniques can be applied to problems with 
cavitating regions. Determining the film reformation 


boundary is more difficult than determining the film 
rupture boundary due to sudden changes in the flow 
variables across the front. 

If the initial location and slope of the boundaries are not 
known, they can be determined by employing a trial and 
error method (the discontinuity fitting method). The loca- 
tion of the boundaries can be assumed and the flow field 
on either side of the discontinuity can be determined. The 
equivalent Rankine-Hugoniot equation can then be ap- 
plied across the boundary to verify the assumption. If the 
initial location and slope are known, then the governing 
equations coupled with the equivalent Rankine-Hugoniot 
equation can be solved to determine the new velocity and 
the new locations of the discontinuities. 

The algorithm introduced by Elrod (198 1) is essentially 
a discontinuity capturing technique. By combining 
the governing equations for the full film and cavitated 
regions and conserving mass flow through the entire 
bearing, the ‘universal equation' is cast into a conservation- 
law form. Hence, the discontinuities can be predicted as a 
part of the solution. This method is simple to implement 
and does not require any knowledge about the location of 
the discontinuities. The boundaries are predicted very 
effectively. 

Concepts from Transonic Flow 
Computation 

The authors have utilized a few transonic flow compu- 
tational concepts in the analysis of cavitated bearings. The 
following is a brief discussion of this work. 

Computational Algorithm 

When the potential flow equation was used in transonic 
flow computation, difficulty was encountered due to the 
reversal of the velocity vector in the supersonic flow 
regime. Murman and Cole (1971), in a landmark paper, 
demonstrated a simple way to obtain a meaningful solu- 
tion by proposing the use of central differencing in the 
subsonic region (elliptic) and one sided upwind 
differencing in the supersonic region (hyperbolic). Jameson 
(1975) created a type differencing scheme by introducing 
an artificial viscosity term into the governing equation. 
This enables automatic switching of the form of 
differencing as required within different regions of the 
flow. Also, with the explicit addition of an artificial 
viscosity term, the conservation form of the equation was 
preserved. 

In the analysis of cavitation in bearings, Elrod (1981) 
modified the originally proposed algorithm by Elrod and 
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Adams (1974), by incorporating the idea of changing the 
form of differencing in the full film and cavitated regions. 
However, this effect was achieved by ‘trial and error* and 
the algorithm was empirically developed. A type 
differencing scheme was developed by Vijayaraghavan 
and Keith (1989), by introducing an artificial viscosity 
term (much like Jameson’s) into the governing equation 
which in turn permitted the algorithm to be mathemati- 
cally derived. In addition, with this modification, the 
discretization was also performed in conservative form. 
The predictions using this modified algorithm were found 
to compare well with the predictions using the Elrod 
algorithm and with experimental data. Hence, the modi- 
fied cavitation algorithm , with this firmer base, was thought 
to offer greater potential for further improvement 

Grid Control 

Numerical grid generation is a fairly common tool to 
model arbitrarily shaped regions in computational fluid 
dynamics. This is basically a procedure to distribute in an 
orderly manner the grid points in the physical field in such 
a way that efficient communication between the points 
and all the physical phenomena on the entire continuous 
field is represented with sufficient accuracy (Thompson, 
1984). Also, the region in the immediate vicinity of solid 
surfaces are dominant in determining the character of the 
flow due the large gradients prevailing in this region. 
Accurate prediction of flow variables in this region is 
important since the final values of the variables strongly 
depends on this boundary prediction. When such high 
gradient regions are not known d priori , a dynamically 
adaptive grid system can be an effective tool. This is an 
active area in grid generation research. By dynamically 
readjusting the grid distribution as the solution proceeds, 
high resolution solutions can be obtained with fewer grid 
points. 

In the case of a cavitated bearing, film rupture and 
reformation locations are not known d priori. Also, accurate 
prediction of the pressure distribution is the primary 
requirement in the full film region, since, all the perform- 
ance parameters depend upon the pressure profile. Hence, 
closely placed grid points around the cavitation bound- 
aries and more grid points in the high pressure gradient 
regions should lead to a more accurate prediction of the 
pressure profile. In the case of a bearing with a misaligned 
journal, the maximum pressure location in the axial direc- 
tion is shifted towards the edge of the bearing and to 
correctly predict the pressure distribution in this region, 
closely spaced grid points are required. With a uniformly 
fine grid arrangement, this will result in a large number of 
grid points located within regions where such a small grid 


spacing does not necessarily provide any significant 
improvement in the accuracy. With a grid adaption tech- 
nique all these effects can be achieved by moving the grid 
points and selectively locating them in such a way that 
accurate solutions can be obtained with fewer grid points. 
The concept of grid generation/transformation and grid 
adaption techniques have been incorporated into the 
modified cavitation algorithm by Vijayaraghavan and 
Keith (1990(a)). 

The grid adaption procedure can be tailored to the 
problem being solved. However, the procedure envisaged 
was to perform initial computations which locate the film 
rupture and reformation fronts, to distribute closely spaced 
grid around them, and then to rearrange the grid in the full 
film region according to the pressure gradient In the case 
of a misaligned journal, when the degree of misalignment 
is large, grid adaption in the axial direction is applied to 
cluster the grid around the maximum pressure location 
(Vijayaraghavan and Keith, 1990(b)). By aligning the grid 
with discontinuities, the flow field is divided into zones of 
full film and cavitated regions. This enables the jump 
conditions to be applied effectively. In addition, in the 
case of a two-dimensional time asymptotic solution, by 
aligning the discontinuity along one coordinate direction, 
according to the equivalent Rankine-Hugoniot equation, 
the mass flow along the other coordinate direction is 
continuous across the discontinuity. This method of grid 
adaption combines features of both shock fitting and 
shock capturing techniques. 

The predicted performance of the bearings using the 
adaptive grid and conventional orthogonal grid arrange- 
ments were found to be comparable. The results obtained 
thus far demonstrate the usefulness of these techniques in 
the analysis of bearing problems. The transformation of 
the governing equation and numerical differencing in a 
nonorthogonal coordinate system could be confidently 
applied, primarily due to the mathematical base provided 
in the modified algorithm. 

Solution Procedure 

In the case of transonic flow, for steady problems, the 
converged solution obtained by using relaxation methods 
generally requires a very large number of iterations due to 
the slow convergence rate. The two most effective solu- 
tion acceleration techniques for rapid convergence in the 
transonic flow computations are approximate factoriza- 
tion of the difference operators and the use of multiple 
grids (Jameson, 1987). For genuine unsteady transonic 
flow problems, Newton iteration techn iques can be applied 
to the governing equation, at every time step, to obtain 
time accurate solutions (Shankar et al, 1985). 
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In the analysis of cavitated bearings, relaxation meth- 
ods are found to also require a large number of iterations 
to obtain a converged steady state solution. Woods and 
Brewe (1989), in the analysis of a dynamically loaded, 
submerged journal bearing, incorporated a multigrid tech- 
nique into the Elrod algorithm and obtained considerable 
savings of computer time. In an unsteady analysis of 
cavitated bearings, particularly dynamically loaded bear- 
ings, time accurate solutions are very important. In such 
cases, the accuracy of the solution can be improved by 
adding a Newton iteration technique. The Newton itera- 
tion scheme and approximate factorization techniques 
were developed by Vijayaraghavan and Keith (1990(c)) 
for the modified algorithm . The approximate factorization 
technique was found to be robust and efficient. The unique 
advantage of these techniques is that with the same unsteady 
formulation, both time accurate unsteady results and fast 
convergent asymptotic steady state solutions can be 
obtained in less computer time than by using a steady state 
formulation. 

Conclusion 

Recognition of the mathematical similarities between 
internal transonic flow and cavitated bearing flow is an 
important and useful concept To the best of authors’ 
knowledge, such similarities have not been pointed out 
before. The analogy has been exploited to determine the 
characteristics within the cavitated region and the jump 
conditions across any discontinuity in the flow field. By 
virtue of similarities between the two flows, advanced 
concepts of transonic flow computations can be incorpo- 
rated into numerical predictions of cavitation in bearings. 
Determination of the reformation boundary using both 
shock fitting and shock capturing methods are discussed. 
With the conservative formulation of the governing ‘uni- 
versal equation’, the shock capturing method is very 
effective and simple to implement. The concepts of tran- 
sonic flow computation have been developed and 
successfully incorporated in three important areas, namely, 
algorithm development, grid arrangement and control, 
and efficient solution computation. The results obtained 
are encouraging and it is believed that many more such 
extensions may be possible which will result in improved 
numerical prediction of cavitation in bearings. 
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